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Abstract 

The plastic response of beryllium was investigated on nanosecond time scales ap- 
propriate for inertial confinement fusion fuel capsules, using laser- induced shock 
waves, with the response probed with surface velocimetry and in-situ x-ray diffrac- 
tion. Results from loading by thermal x-rays (hohlraum) were consistent with more 
extensive studies using laser ablation. Strong elastic waves were observed, up to ~1 
km/s in free surface speed, with significant structure before the arrival of the plastic 
shock. The magnitude and shape of the precursor could be reproduced with a plas- 
ticity model based on dislocation dynamics. Changes in lattice spacing measured 
from the x-ray diffraction pattern gave a direct measurement of uniaxial compres- 
sion in the elastic wave, triaxial flow from the decay of the precursor, and triaxial 
compression in the plastic shock; these were consistent with the velocity data. The 
dynamic strength behavior deduced from the laser experiments was used to help 
interpret surface velocity data around the onset of shock-induced melting. A model 
of heterogeneous mixtures is being extended to treat anisotropic components, and 
spall. 
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1 Introduction 



The dynamic plasticity of beryllium on nanosecond time scales is important 
in the development of inertial confinement fusion (ICF) ([21 . The ini- 
tial experiments on thermonuclear ignition at the National Ignition Facility 
will use 192 laser beams to deliver 1-2 MJ of energy inside a hohlraum. The 
resulting flux of soft x-rays will implode a spherical capsule containing cryo- 
genic deuterium-tritium (DT) fuel. If the implosion is sufficiently vigorous and 
uniform then thermonuclear burn is predicted to propagate from the central 
hotspot. The implosion has a relatively high aspect ratio (initial to final radius) 
of around 30, and is unstable with respect to spatial variations in thickness, 
surface finish, or radiation drive. 

A leading material for the fuel capsule is Be containing up to around 1 weight 
percent of Cu to increase its opacity to the radiation drive. Compared to al- 
ternative capsule materials such as plastic, Be is dense and has a high thermal 
conductivity, which helps maintain a uniform layer of DT ice. However, Be ex- 
hibits elastic and plastic anisotropy when subjected respect to mechanical or 
thermal loading. This anisotropic response is likely to induce velocity fluctua- 
tions from the hohlraum drive, which may seed significant implosion instabili- 
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ties. It is thus necessary to define a specification for the microstructure - e.g. a 
maximum grain size - for Be-Cu in an ICF capsule. The response of Be on rel- 
evant length and time scales, ~100 fim and ~10ns, has not been characterized 
previously, so this work includes the development and calibration of appro- 
priate models, continuum mechanics simulations in which the microstructure 
is resolved in order to predict the velocity fluctuations directly, and the ex- 
perimental validation of predictions of anisotropic response and instability 
seeding. The continuum mechanics simulations are computationally expensive 
in three dimensions, so some systematic studies of the velocity fluctuations 
are performed in fewer dimensions. A mixture model, originally developed for 
reactive flow studies (J9[) but subsequently used for other applications including 
the shock compression of porous minerals (J6|), was generalized to predict the 
average response of strong, polycrystalline materials, as an inline treatment 
of texture at coarser resolutions. 

This paper discusses experiments performed to investigate the elastic-plastic 
response of Be on nanosecond time scales, and the development of physically- 
based models to describe the behavior. The mixture model, with extensions 
for strength, is presented. The treatment of strength should allow material 
damage and failure to be modeled seamlessly; damage pertinent to dynamic 
failure by spallation is discussed. 



2 Ablative loading experiments 

Plastic flow in Be under dynamic loading was investigated experimentally 
using laser ablation to induce shock waves. Experiments were performed at the 
TRIDENT laser facility at Los Alamos, which has been used to induce shocks 
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and quasi-isentropic compression waves by ablation of a free surface ( 11 71 ; Il5l ) 



and also to induce loading and accelerate flyer plates by ablation through a 



transparent substrate to act as confinement (Il6l ). TRIDENT has a Nd:glass 
laser system with a fundamental wavelength of 1054 nm. For more efficient 
coupling with matter, the pulse was frequency-doubled to 527 nm (green). 
Operating in nanosecond mode, the pulse comprised up to 13 elements 180 ps 
long for each of which the intensity could be set independently, allowing the 
irradiance history to be controlled, with a total energy of up to 250 J. In each 
experiment, the planar sample was irradiated on one side, using a Fresnel zone 
plate was distribute the laser energy uniformly over a spot 5 mm in diameter. 
(Fig. 2.) 

Velocity histories were obtained at the surface opposite the laser ablation, 
by laser Doppler velocimetry over a line 600 /im long centered opposite the 
center of the ablation spot. In some experiments, a second laser beam was 
used to generate x-rays by focusing it over a ~150 /im spot on a Ti foil, 10 fim 
thick; the plasma resulting from a pulse energy of 100-200 J in 1-2.5 ns cooled 
by emitting predominantly He-a-like line radiation (Fig. 3). The x-rays were 
used to obtain diffraction measurements in situ during the shock wave event 



and thus to monitor the response of the crystal lattice directly (Il2l ) 



To simplify the analysis and interpretation of shock wave experiments, it is 
desirable that the loading history should be constant for some period, then 
fall rapidly to zero. A disadvantage in using ablative loading compared with 
flyer impact experiments is that the relation between irradiance history and 
pressure history applied to the sample is not simple: for instance, a constant 
irradiance does not lead to a constant drive pressure jl7). The irradiance his- 
tory must in general be tailored to the material, desired pressure, and pulse 



4 



duration. Often on nanosecond time scales, a constant pressure requires an 
irradiance history which increases by several tens of percent over the duration 
of the pulse. Tailored pulse shapes such as this can be generated at TRI- 
DENT because of the flexibility in pulse construction, but would generally 
require several iterations because of inaccuracies in radiation hydrodynamic 
simulations used to predict the desired shape, making it necessary to per- 
form experiments on shock flatness to verify and adjust the shape. The pulse 
shape and energy produced by the laser were not perfectly reproducible, so 
there is a limit to the accuracy with which the pressure history can be spec- 
ified - e.g., constant - on a given shot. Rather than attempting to use an 
ideal irradiance history, the pulse shape was measured on each shot and used 
in radiation hydrodynamics simulations to predict the pressure history at a 
point close to the ablated surface. The pressure history was then used as a 
time-dependent boundary condition in continuum mechanics simulations, al- 
lowing more detailed models of the behavior of condensed matter to be used 
without requiring the simulation program to include radiation physics also. 

ications on 
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This procedure has been found to work well for a range of app 
metallic and non-metallic elements, and intermetallic compounds (117 
largely because relevant material properties in the plasma state - the equation 
of state (EOS), opacity, and conductivities - can be predicted with adequate 
accuracy. For the pressures of a few tens of gigapascals considered here, the 
radiation hydrodynamics simulations predicted that less than a micrometer 
of material was ablated. Recovered samples have shown only a thin (~1 /im) 
layer of the remaining material significantly affected by the heat associated 
with the ablation process. 



At a few tens of micrometers, the sample thickness was two orders of magni- 
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tude less than the spot size, so the experiments were one-dimensional during 
the passage of the initial shock and subsequent release from the drive surface 
and the surface monitored by velocimetry. One-dimensional simulations of 
radiation hydrodynamics and continuum mechanics were adequate. The radi- 



ation hydrodynamics simulations used the computer programs 'LASNEX' ( I22l ) 
and 'HYADES' (jj). The continuum mechanics simulations used the program 
'LAGC1D' (1201 ) . These programs used second-order finite difference represen- 
tations of space and time with Lagrangian (material-following) derivatives, 
and artificial viscosity to stabilize shocks. 



3 Shock response data for Be 

The plasticity experiments on Be included single crystals and rolled foils. The 
crystals were cut nominally parallel to (0001) planes from a boule grown by 
zone refining, and were polished to 30-40 fim thick. The rocking curves were 
around 2° wide, full width, half maximum (Fig. 4). The rocking curve was 
determined statically by diffraction of Cu K a x-rays, integrated over almost 
the whole of each sample. The profile of the dynamic diffraction lines was 
consistent with the static rocking curve, though there was some sign that in 
the dynamic diffraction geometry, where the source was quite close to the 
sample, the width of the diffracted lines was slightly narrower. The foils were 
commercial high-purity material (at least 99.8% Be). The texture, measured by 
x-ray orientation mapping, was found to be with c-axes oriented predominantly 
normal to the foil - i.e., similar to the single crystals but with a distribution 
of angles ~ 20° wide (full width) - and the a-axes distributed continuously in 
the plane, though clearly exhibiting preferential directions (Fig. 5). Compared 
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with the single crystals, the foils presumably had been subjected to a much 
greater degree of plastic work. 

In the velocity history records from the single crystal experiments, the elas- 
tic precursor was very clear with the lowest pressure drive (mean pressure 
15.0 GPa). The velocity showed a sharp initial peak (rise time faster than 
could be resolved), followed by a deceleration, followed by the main plastic 
wave (Fig. 6). The x-ray diffraction records for the three lowest pressures ex- 
hibited a peak corresponding to the uniaxial compression in the precursor 
deduced from the initial peak, and at the lowest two pressures an additional 
peak corresponding to the same volumetric compression but expressed isotrop- 
ically rather than uniaxially - suggesting that plastic relaxation to an isotropic 
state occurred before the plastic wave arrived (Figs 7 and 8). Velocity histo- 
ries from experiments on rolled foils of similar thickness showed a precursor 
of lower amplitude, with a finite and reproducible rise time of around 400 ps 
(Fig. 6). The amplitude of the precursor varied more between experiments, 
with higher amplitudes observed in thinner foils. The higher amplitudes could 
be caused by time-dependence of plastic flow, or by greater plastic strain im- 
parted by rolling to a thinner foil. 

The laser Doppler diagnostic recorded the velocity history over a line across 
each sample. The records obtained from single crystals were spatially uniform, 
demonstrating that the laser drive was adequately smooth over the surface of 
the sample. Some of the foil records exhibited spatial variations in velocity 
following the precursor, in particular regions in which there was a deceleration 
after the initial elastic peak, resembling the single crystal velocity history 
(Fig. 9). 
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These velocity data - particularly the single crystal record - were compared 
against simulations using different models of dynamic plasticity. There was 
some uncertainty in sample thickness, and there was also some timing jitter 
in the signal used to trigger the diagnostic cameras, so timing variations of up 
to ~ 1 ns were plausible. 



4 Simulations with different plasticity models 



Several models have been developed to describe the plastic response of Be 



under dynamic loading (11 ll ). based on gas gun experiments using samples 
of the order of millimeters thick. The laser-driven samples were a few tens 
of micrometers thick, so any time-dependence in the plastic flow processes 
would be expected to lead to a higher flow stress and thus a higher velocity in 
the elastic precursor; this is in fact what was observed. One of the published 
models (Steinberg-Lund) is time-dependent, but these experiments explored 
time scales far shorter than the regime in which the model was calibrated, so 
one would not expect it necessarily to perform well. 

The speed of the free surface from the elastic precursor was a significant frac- 
tion of the speed induced by the plastic shock. The elastic wave reaches the 
surface first, reflects from it, and rattles between the surface and the shock. 
For a given shock pressure, the maximum speed seen at the free surface de- 
pends on the amplitude of the elastic wave. For a perfect model of the drive 
pressure, one would not expect the peak free surface velocity to match the 
experiment unless the elastic wave was modeled adequately. 

Spall was ignored in most of the simulations, so the calculated surface velocity 
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decelerated close to zero after the end of the drive pulse. 



4-1 Hydrodynamic shock 



For reference, simulations were performed with no strength model, i.e. the 
scalar EOS only. A cubic Griineisen form was used (111: it was found previously 
that all common EOS agree for the mechanical response of Be in the regime 
explored by these experiments. With no elastic wave, the peak surface velocity 
was higher than observed; the arrival time of the shock was consistent with 
the experiment. (Fig. 10.) 



4-2 Elastic-perfectly plastic 



As the plastic behavior of Be had not previously been explored in the time 
and length scales investigated here, it is prudent to perform simulations using 
a model with as few parameters as possible. In the elastic-perfectly plastic 
strength model, the material state is described by the isotropic state and the 
elastic strain deviator. The constitutive model comprises the EOS, a shear 
modulus G, and a flow (or yield) stress Y . The elastic response is a = Ge e ; if 
the norm of the stress ||cr|| < Y then e p = 0; otherwise e p = e a . 

Parameters reported for beryllium are G = 151 GPa, and Y between 0.33 and 
1.3 GPa. With values in this range, the amplitude of the precursor was much 
lower than observed experimentally. The initial peak could be reproduced 
with Y ~ 6.0 GPa and the following trough with 4.0 GPa; the shape was 
not reproduced using this model. With high values of the flow stress, the 
effect of the reverberating elastic wave was evident as a sequence of steps in 
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the acceleration history. The steps were similar to some features seen in the 
experimental velocity history. As the flow stress was increased, the predicted 
peak velocity decreased to approach the measured value, and the release time 
increased to be close to the observation. (Fig. 11.) 



4-3 Steinberg- Guinan and Steinberg-Lund 



The Steinberg-Guinan model (jlll ) is a generalization of the elastic-plastic 
model to include work-hardening and thermal softening, and varying the ini- 
tial plastic strain e p o to investigate the sensitivity. No set of parameters was 
found which was capable of reproducing the detailed structure of the precur- 
sor. If the maximum flow stress was increased to give an approximate match 
to the amplitude of the precursor, the amplitude and duration of the plastic 
peak became very different to the observed velocity record. 



The Steinberg-Lund model (jlll ) is a generalization of the Steinberg-Guinan 
model to include time-dependence. Using the published parameters for Be, the 
amplitude of the elastic wave was lower than observed and the plastic shock 
released much earlier then observed. No set of parameters was found which 
reproduced the observed velocity history. 



4-4 Simplified micro structural plasticity 



Rather than attempting to test and calibrate explicit microstructural plasticity 
models directly against the experimental data, simple extensions of the elastic- 
perfectly plastic model were made to include time dependence with as few 
other parameters as possible. Thus the systematic behavior observed in the 
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single crystal experiments could be investigated in isolation from the effects 
of different terms in the more complicated models. 

The plastic strain rate was chosen initially to be 



where the shear modulus G, activation temperature T*, and rate multiplier 
a are in principle functions of the mass density p, but were kept constant for 
the low compressions studied here. The activation temperature was initially 
chosen to match the Steinberg-Lund value, T* = 3600 K, and the sensitivity to 
a and T* was explored. It was possible to reproduce the observed wave profile 
qualitatively; in particular, deceleration could be induced following the ini- 
tial elastic peak. However, it was not possible to reproduce the precursor and 
the rise of the plastic wave very accurately. The simulations exhibited many 
small-amplitude spikes, caused by under-damped reverberations of the pre- 
cursor. Work-hardening from plastic flow behind the initial wave would cause 
a smearing of the plastic wave and should act to damp out the reverberations, 
but this was not investigated in the context of this model. 

A slightly more realistic flow model was investigated to add the flexibility to 
reproduce the characteristic decay time of the precursor independently of its 
initial amplitude. In reality, the barrier for defect hopping mechanisms in a 
crystal is altered by the strain. The potential of an atom or defect with respect 
to displacement is periodic, with an amplitude like sin 2 (7re/6) where e is the 
linear displacement and b the length of the Burgers vector in the direction 
of interest. Thermally-activated hopping may in principle occur in either the 
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direction of increasing or decreasing strain energy; the net rate is 
e p = a{p) {exp[-T*(l - 7 )/T] - exp[-T*(l + 7 )/T]} 



(2) 



where 

-*.) - sin* (^) ■ (3) 

The reference period for elastic strain was taken to be a constant, b(p) = 1. 
The activation temperature was again chosen to match the Steinberg-Lund 
value, T* = 3600 K, and the sensitivity to a was explored. For a < 10/ns, the 
material propagated an elastic wave. At a = 100/ns, the shock had started 
to split into a precursor and a plastic wave, though elastic reverberations in 
the plastic wave were pronounced. For a > 1000/ns, the precursor had the 
same peak-and-dip shape as seen experimentally, with a step just less than 
halfway up the plastic wave. A good fit to the initial peak and the subsequent 
deceleration was obtained for a = 2000/ns, though the acceleration during the 
plastic wave was sharper. As with the plastic relaxation model, overall rise time 
of the plastic wave matched the experiment because of the delay introduced 
by the step during the rise. The release wave after the peak occurred at the 
same time, and some structures in the peak resembled the experiment. 

For defect-induced plastic strain, the rate is 

e = R d \b\ 2 l 

where I is the length of dislocations per volume of material, |6| the Burgers 
vector (projected onto the direction of interest), and the jump rate of a 
defect. One could assume that a defect jog must comprise N atoms in order 
to be stable, where N could be calculated given the line energy of the defect 
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compared with the strain energy reduced by motion. If the probability rate 
of a single atom jumping is R a , the probability rate of a stable jog forming is 
Rq . As a simple estimate, the local jump rate R a occurs over a length |5|, so 

e = R a \b\L 

For Be, |6| ~ 2 A. The dislocation density can be estimated from the rocking 
curve, width 6. Assuming the curvature is caused by a population of defects of 
the same orientation, the linear spacing between defects is \b\/8, so the defect 
density 

m 

For a 2° rocking curve in the Be crystals considered here, this implies I ~ 
1.7 x 10 10 /m. The rate at which individual atoms attempt to jump and thus 
mediate motion of a defect is ~ 10 12 Hz, thus the prefactor a = 2000/ns is 
consistent with the rate of deformation that one would expect. 

Spallation of the sample should really be included through a failure model 
intimately connected with the plasticity model. Here, simple spall simulations 
were performed to test whether the results would be at all consistent with 
the experiment, given the strong and structured elastic precursor. Spall was 
modeled using the 'minimum pressure' method: the EOS was modified so that 
the minimum pressure it could produce was limited to a specified value - 
the spall strength, defined as a negative number - and any greater degree of 
tension would cause unbounded expansion, leading to a region of low density. 
The degree of deceleration after the plastic peak was reproduced with a spall 
stress ~2 GPa, compared with a reported value from gas gun experiments of 



0.35 GPa for Be alloy S200 (jlll ). The difference is reasonable on grounds of 
rate-dependence. (Fig. 12.) 
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4-5 Micro structural plasticity 



The simple models discussed above are adequate to demonstrate that mi- 
crostructural processes explain the plastic behavior of a single orientation of 
Be crystals, but a more complete model of rate-dependent crystal plasticity is 
necessary for the resolved-microstructure simulations of instability seedin g fo r 
ICF; the discussion below follows the spirit of Rashid and Nemat-Nasser fl 1 Of ) . 
The crystal lattice is described locally by the lattice vectors, {a,b,c}; these 
can be defined with respect to standard conventions of orientation (a parallel 
to the x-axis, b in the xy plane) by a local rotation matrix, R. The EOS and 
elasticity give the stress tensor s(a,b,c). Plastic flow is mediated by defects 
- dislocations or disclinations - which act along a set of slip systems defined 
with respect to the lattice vectors. 

Each slip system is characterized by a plane (normal n) and a direction d in 
the plane. Given the stress tensor s, the stress resolved along this slip system 
is 

s r = d.s.n = d.a.n (4) 
since n and d are orthogonal. The plastic flow rate in the system is 



e r = Z(p, T) {exp [-T*(l - 7 )/T] - exp [-T*(l + T )/T]} f(p d / Po ) (5) 



where pd is the dislocation density, p Q the obstacle density, T* the barrier 
energy (again expressed temperature), and 




(6) 
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where 



d.e.n = d.e.n. (7) 



Given a plastic strain rate e r in a system, the strain tensor changes as 
e = e = — V] e r d ® n 



(llOl ). Note that e is always symmetric, though contributions from individual 
slip systems may not be. So far, the strain rate in a given slip system is identical 
with the simple hopping-based model above except for 0(1) geometrical terms 
describing the resolved shear stress on the system and the resolved strain, and 
the explicit magnitude in terms of density of dislocation length. Further elastic 
precursor measurements on Be crystals of different orientations will be needed 
for a complete calibration of the crystal plasticity model. 



Work-hardening enters crystal plasticity through the interactions between dis- 
locations or twin boundaries and obstacles, including other dislocations. The 
dislocation density and obstacle density increase with plastic strain: 

Pd = F l {p,p )e r + F 2 (cr r ), p = G(p,p )e r (9) 



where F\ and F2 are dislocation generation by flow (e.g. Frank-Reade sources) 
and by direct nucleation respectively. These parameters will be calibrated 
against molecular dynamics simulations, and validated by comparison with 
the slope and detailed shape on the main plastic wave. 
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5 Strength of heterogeneous mixtures 



A large amount of work has been done on interpreting shock wave data on 
heterogeneous materials - particularly porous ones - in terms of properties 



of pure materials (Il9l ). Work has also been done to predict restricted sets of 



mechanical properties such as the shock Hugoniot, given the composition or 



a. 



initial porosity (Il9t y). However, relatively little has been done to develop 
models capable of predicting wide-ranging response in a form satisfying the 
requirements of general continuum mechanics simulations. 

It is highly desirable to have a model which predicts the response of a het- 
erogeneous material given models for homogeneous components and a repre- 
sentation of the microstructure. For mixtures and materials of low porosity, 
models of the response for use in continuum simulations have been developed 
from experiments on the heterogeneous material itself, without including a 

(1; 18; 7) 



model of the "pure" components 



7 ). This approach is pragmatic, but 



has significant disadvantages. If interest in a different composition grows, fur- 
ther experiments are needed to calibrate a new model. Ab initio techniques are 
increasingly capable of predicting the response of pure elements and of stoi- 
chiometric compounds; these theoretical models cannot be applied to mixtures 
and porous materials without a predictive model of heterogeneous materials. 
Finally, for some applications, the composition or microstructure of a hetero- 
geneous material may vary continuously in space, or may evolve in time. This 
is the case for sputtered Be-Cu alloy in ICF, where the texture may change 
radially through the capsule, particularly in the presence of an ion current. 
A predictive, inline mixture model is desirable here for simulations at varying 
degrees of detail from resolved microstructure to continuum-averaged proper- 
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ties. 



We have previously developed explicit mixture models for chemical explosives, 
allowing the shock initiation properties to be predicted as a function of compo- 
sition and porosity (13; 0). These models were essentially hydrodynamic, i.e. 
they did not incorporate elastic and plastic effects explicitly in the mixture 
model or the stress field generated by the material. Plasticity was included 
as a contribution to heating in the treatment of hotspots. Here we present a 
more general model in which the constitutive behavior of the components is 
also included. 

The dynamic response of a heterogeneous mixture is determined by the prop- 
erties of its components and the equilibration of stress and temperature be- 
tween them. For the propagation of a single shock wave of constant pressure, 
the conditions of interest may be perfect equilibration. However, there is no 
such thing as an 'ideal' shock wave experiment: in practice, experiments and 
real applications explore shock propagation through a finite thickness of ma- 
terial and hence probe a finite time scale. If a steady shock is shown to exist, 
then the Rankine-Hugoniot equations apply and equilibrium conditions ap- 
ply behind the shock. In practice, the steadiness is established with a finite 
accuracy, and in general there may be equilibration processes occurring on a 
longer time scale which are in effect frozen during the shock wave experiment. 
Thus it is not always correct to assume that perfect equilibration occurs across 
a shock wave. In the case of more general loading histories, the response of 
the material may depend on its instantaneous non-equilibrium state, so this 
should be taken into account. It has been found by comparison with shock 
experiments on mixtures that mixture models which ignore equilibration and 
strength effects are not able to reproduce all the observed shock behavior 
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(1l9l). Thermally-activated processes such as chemical reactions, plastic flow, 
and phase transitions depend more strongly on non-equilibrium effects. 



The model described below was developed to act as a bridge between simu- 
lations with a resolved microstructure and continuum-averaged simulations. 
It is based on continuum- level averaging, but includes microstructural terms 
which can be estimated fairly easily, or calibrated against more detailed sim- 
ulations. Each component is described as a pure material, i.e. with its own 
EOS and strength model. The model was designed for use in continuum me- 
chanics, so it consists of a local state (varying with location and time inside 
a material subjected to dynamic loading) and a description of the material 
response (the same throughout a particular material). Since each component 
has its own state, a stress and temperature can be calculated for each. The 
material response model uses the local state to calculate 'external' properties 
such as mass density, stress and temperature, and to calculate the evolution 
of the state given applied loading and heating conditions. Stress and tempera- 
ture are equilibrated explicitly according to a separate time constant for each. 
This produces an exponential approach to mechanical and thermal equilibrium 
separately. Ideal equilibration can be enforced by setting the time constants 
to a small value. Stress equilibration was performed by adjusting the strain 
state of each component towards equilibrium. In the absence of constitutive 
effects, stress equilibration becomes pressure equilibration; this was performed 
by adjusting the volume fraction of each component towards equilibrium, and 
expanding or compressing along an isentrope. Temperature equilibration was 
performed by transferring heat energy at constant volume. Equilibration based 
on finite rates was found to be convenient because it can be designed to use 
essentially the same material properties as a continuum simulation of a 'pure' 
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material, thus is was possible to implement the heterogeneous mixture model 
as a superset of existing types of material property (EOS etc) without having 
to add additional types of calculation (e.g. different derivatives of the EOS). 
The numerical schemes are described in greater detail below. The time scales 
for equilibration may be estimated from the grain size and properties of each 
component. In principle, they could be considered as additional continuum 
variables and evolved according to other microstructural processes. (Fig. 13.) 

Porous materials were represented by starting with a non-zero volume fraction 
of gas (air or reaction products). 



5. 1 Representation of the micro structure 



We have previously reported results from reactive flow models using mixtures 
of EOS (jl3l ; [8|; |9j). These models treated the mixture as a set of discrete 
components i of volume fraction f i} each with its own state Sj. This approach 
can be extended in a straightforward way to equilibrate stresses rather than 
pressures. The representation of the microstructure was contained entirely in 
the volume fractions. Here we describe an incremental improvement describing 
the microstructure with slightly more detail. 

For each component, a 'fractional overlap area' uji was also defined; this tensor 
quantity expresses the connectivity of the component across the mixture in 
each direction, in the same sense as stress and strain tensors. The overlap 
areas were used in equilibrating stress, allowing a strong component deemed 
to be continuous across the mixture to support a stress in the presence of 
another component - such as a void - which could accommodate strain with 
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a smaller stress (Fig. 14). An overlap area of zero for a component in some 
orientation implies that the component is continuous across the mixture in 
that orientation. An overlap area of unity implies that the component stress 
should be equilibrated ideally in that orientation. Intermediate values were 
used to simulate granular mixtures. 

Each component in general has different material properties. In addition to 
the properties pertaining to each 'pure' component, each component was asso- 
ciated with a tensor evolution parameter dui/de describing the rate of change 
of the overlap area with plastic strain normal to the overlap orientation. This 
parameter was used typically to simulate the homogenization of the mixture as 
shear strain breaks the continuity of components across the mixture. The evo- 
lution of the overlap areas was approximated by a material constant (vector) 
dujjdt where e is interpreted as the strain component normal to the compo- 
nent of uj. A stress tensor and temperature is calculated for each component, 
and for the mixture as a whole. 

5.2 Thermal equilibration 

Thermal equilibration acts to reduce differences between the temperature 7* in 
each component, which is a function of the instantaneous state Sj in that com- 
ponent. Given the component temperatures, we define a mean temperature 
for the mixture as a whole, weighting by the heat capacity of each component: 



where /j is the volume fraction of the ith component, c v is its specific heat 
capacity, and p its mass density. Given the thermal relaxation time scale tt, 




(10) 
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the relaxation factor over a time increment At is 



4>r= (l - e' At ^) . 



(11) 



Thus the change in specific internal energy for each component can be calcu- 
lated for a finite time increment, 



5.3 Mechanical equilibration 

Mechanical equilibration acts to reduce differences between the stress Tj in 
each component, which is a function of the instantaneous state s« in that com- 
ponent, which includes the strain. Here, stress and strain are tensor quantities, 
expressed in Voigt (vector) notation if desired. Given the component stresses, 
we define a mean stress for the mixture as a whole: 



Given the mechanical relaxation time scale r m , the relaxation factor over a 
time increment At is 



In each component, the difference in stress from the mean can be expressed 
as a strain increment Atj toward equilibrium, 






(13) 




(14) 




(15) 
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where S is the stiffness matrix for the component's material in its instanta- 
neous state, 




d[Ti] k 



(16) 



d[e]j ' 



The strain increment can be expressed as a strain rate, 



gradw 



(17) 



At 



and the change in volume fractions is 



dfi = l 
dt 3 




) 



(18) 



In the case of a material described by a scalar EOS, the equilibration relations 
can be expressed in simpler form. Now equilibration acts to reduce differences 
in the component pressures Pi, functions of the instantaneous states Sj. In each 
component, the difference in pressure from the mean can be expressed as a 
volume change A/j toward equilibrium, 



where the mean pressure p is the volume-weighted average of the component 
pressures. The isotropic stiffness is 




(19) 



dpi = Cj(sj)pi 
df* ~ fi 



(20) 



where q is the bulk sound speed. 
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6 Spall as extension to heterogeneous strength 



The heterogeneous strength model described above is a natural framework for 
developing a class of models for material failure, including spall. Key elements 
are the inclusion of void as a component in the mixture, and the use of the 
overlap tensor to represent the orientation of growing regions of void with 
respect to the direction of maximum tension. The overlap tensor allows the 
model to discriminate between ductile and brittle failure (Fig. 15). 



The formation and propagation of cracks and voids is dominated by stress 
concentrations, i.e. localization phenomena caused by the heterogeneity of the 
microstructure. Stress concentration is not treated adequately by the overlap 
tensor, so sub-scale physics is required, at least at the level of an effective flow 
stress for the solid microstructure under tension. This is similar to the situation 
with heating induced by pore collapse in the previous work on heterogeneous 
high explosives 9), where additional shape functions were used to model the 
enhanced heating around a collapsing pore, i.e. representing the width of the 
temperature distribution as well as the mean. The heterogeneous model should 
not be expected to predict the explicit formation and evolution of pre-existing 
or incipient damage sites. 



The model seamlessly treats recompaction of a damaged or spalled layer, but 
may exhibit unphysical healing, unless the model of texture evolution is so- 
phisticated enough. 
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7 Conclusions 



Elastic-plastic data were obtained for Be (0001) crystals using laser-driven 
ablative loading, with diagnostics comprising laser Doppler velocimetry and 
in-situ x-ray diffraction. These experiments explored the response of Be to 
dynamic loading on nanosecond scales. The flow stress was much higher than 
has been observed on microsecond scales, and exhibited relatively complicated 
temporal structure. Important aspects of the velocity history were reproduced 
using a model of plastic flow which was based on microstructural processes, 
forming the basis of a complete model of crystal plasticity in Be. A hetero- 
geneous mixture model developed previously for reactive flow studies was ex- 
tended to treat strength with a representation of material texture, and seems 
to be a promising platform to treat spall. Accurate models of plasticity and 
phase changes in the homogeneous components of the mixture are needed to 
support physics-based damage and spall models. Such physics-based models 
are necessary to relate the large pool of experimental data and theoretical 
work from microsecond time scales and slower to the qualitatively and quan- 
titatively different behavior observed in high energy density experiments. 
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Fig. 1. Frame from illustrative simulation of shock propagating through resolved mi- 
crostructure in Be. Upper figure shows grains (colored according to the crystal orien- 
tation) and the computational mesh. Lower figure shows thermal colormap of mean 
pressure (black: zero, white: 25GPa). The computational domain is 100 x 50 /im, 
driven by a constant 20 GPa from the lower surface. An elastic precursor (red) can 
be seen running ahead of the main shock, and spatial variations are evident across 
the width. Velocity spectra can be extracted from these simulations and used to 
predict the seeding of implosion instabilities. 
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Fig. 2. Schematic cross-section of laser drive experiments. The sample was 10 to 
50 fj,m thick; the focal spot was 5 mm in diameter. 
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Fig. 3. Schematic of experimental configuration for transient x-ray diffraction mea- 
surements. The x-ray source was 10-15 mm from the sample, and the snouts of the 
x-ray detectors were 50 mm apart. 
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Fig. 4. Rocking curves (distribution of plane orientation with respect to surface) for 
nominally (0001) beryllium crystals. Numbers are TRIDENT shot indices. 
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max. =10. 80 




Fig. 5. Orientation imaging maps of rolled beryllium foil. Each circle shows the dis- 
tribution of orientation of a particular Bragg reflection as a stereographic projection 
centered on the normal to the foil. The color map shows the density of the reflection 
at each angle, on a logarithmic scale. The plotting software used primitive plane 
indices, so '100' refers to the hexagonal plane (1010) etc. 
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Fig. 6. Free surface velocity history from single crystal (TRIDENT shot 12198: 
40 Aim (0001), 2.7PW/m 2 for 1.8 ns) and foils (shot 12182: 50jtm, 3.2PW/m 2 for 
1.8ns; shot 12184: 50 fim, 6.5PW/m 2 for 1.8ns). Foil histories were displaced in 
time so the precursors occurred at a similar time to that in the crystal. 
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Fig. 7. Time-integrated x-ray diffraction data from single crystal and fit by de- 
convolving using measured rocking curve (TRIDENT shot 12202: 40 /um (0001), 
6PW/m 2 for 1.8 ns). Lines were identified by inspection and by comparison with 
compression predicted using radiation hydrodynamics simulations. 
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Fig. 8. Time-integrated x-ray diffraction data from single crystal, converted to vari- 
ations in lattice parameter (TRIDENT shot 12202: 40 /mi (0001), 6PW/m 2 for 
1.8 ns). Width of each bar indicates uncertainty. 
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Fig. 9. Free surface velocity history from a beryllium foil, with traces taken at differ- 
ent points across the sample - field of view ~ 600 /iin - showing spatial variations 
in the precursor (TRIDENT shot 12145: 13 /*m foil, 0.6PW/m 2 for 3.6ns). This 
experiment was driven at a relatively low pressure (~ 7GPa), so the plastic wave 
was only slightly larger than the precursor; the precursor was also stronger than in 
the other foils, presumably because of greater cold working and/or time-dependence 
of flow. 
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Fig. 10. Comparison between free surface velocity history measured for ablative 
loading of a single crystal cut parallel to (0001) (TRIDENT shot 12198: 40 fim 
(0001), 2.7PW/m 2 for 1.8ns) and a hydrodynamic simulation in which spall was 
ignored. 
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Fig. 11. Comparison between free surface velocity history measured for ablative 
loading of a single crystal cut parallel to (0001) (TRIDENT shot 12198: 40 ^m 
(0001), 2.7 PW/m 2 for 1.8 ns) and elastic-perfectly plastic simulations using different 
values of the flow stress Y. 
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Fig. 12. Comparison between free surface velocity history measured for ablative 
loading of a single crystal cut parallel to (0001) (TRIDENT shot 12198: 40 fim 
(0001), 2.7PW/m 2 for 1.8ns) and simulations using the crystal relaxation model, 
(T* = 3600 K and a = 2000/ns), with a simple spall model, for different values of 
the spall strength. 
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Fig. 13. Schematic of the heterogeneous model. 
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Fig. 14. Schematic of the area overlap model. Outer rectangle represents a spatial 
domain described by the mixture model, e.g. a hydrocode cell. Different components 
of the mixture are represented by rectangles of different color. Overlap tensor u 
describes the extent to which each component spans the whole domain, in a given 
direction. 





Fig. 15. Schematic of different types of void growth as represented by the overlap 
model, under strain in the horizontal direction: ductile (left) and brittle (right). 
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Thermal: time scale 

w-T- 

Mean temperature T = — — : wi = fic v {si)p{si). 

Relaxation factor over time increment At: <j> r = (l — < 
Change in specific internal energy: Ae^ = —4> r c v {si) (1 



Mechanical: time scale r m 
Mean stress: f = /^(sj 

Relaxation factor over time increment At: 4> r = (l — . 

Strain increment: Ae^ = r ^(r^ — f) 

where 5 is stiffness matrix: [S^j = dfo^/dlXlj. 

Thus strain rate gradix = Ae^/At, 

and volume fraction change dfa/dt = iTr (dei/dt) . 



